% Process the eej strength derived from satellite data
% and compare it with ground deribed eej indices
%start date 30 OCT 2005
%latest date 21 JUNE 2006
%The latest changes were done for Stefan's new set of data sent on 20-6-6
clear;
load c:\manoj\projects\eej\EEJ_Stefan data;

load c:\manoj\projects\eej\OBSdH;
%OBSdH1 description 
%The file is 15x1096x24 size (15 stn 1096 days 24 hours - dH)
%[1-ABG,2-ETT,3-HYB,4-TIR,5-HUA,6-FUQ,7-MBO,8-GUI,9-AAE,10-QSB,11-GUA,
%12-CBI,13-PND,14-ELT,15-AAE(new),16-NGP]
load c:\manoj\projects\eej\mjdkpf107 mjd kp;
%-------- correlate the CHAMP derived EEJ strength and ground delta H
% with LT correction
%--------------------------------------------------------------------

bin5 = [-60 -55 -50 -45 -40 -35 -30 -25 -20 -15 -10 -5 0 5 10 15 20 25 30 35 40 45 50 55 60];
bin10 = [-60 -50 -40 -30 -20 -10 0 10 20 30 40 50 60]; 
bin20 = [-60 -40 -20 0 20 40 60];

%Control Center

obsdata = 11;  %ETT-HYB = 1,ETT-ABG = 2;TIR-HYB = 3;TIR-ABG = 4; HUA-FUQ 5; MBO-GUI = 6; 
%AAE-QSB 7; GUA-CBI 8; YAP-BIK 9; PND-HYB 10; ETT-PND 11; TIR - PND 12;
%AAE-ELT 13
satdata = 2; % %max current = 1; av cur distan 2;
obslt = [7.9 16.1]; %NOTE the LT are represented as for eg 8.99999 or 9.001 the filter
incrlat = 10;
x1data = bin10;
mindate = tmjd(2000,01,01,0,0,0,1); % starting date
maxdate = tmjd(2002,12,31,0,0,0,1)+1; % ending date
sqcor = 0; %SQ correction
ltcor = 1; %LT correction
robft = 1; %Robust outlier elimination (just remove the most extreme datum)
plotopt = 'b-';
plotcol = 'b';
plotchar = 'b';
ncoeff = 9;%same for hourly means 20 for minute means
kp_thrashold = 2.0;

%load c:\manoj\projects\eej\guacm4.mat index sq_est

switch obsdata,
    case 1,%ETT-HYB
        %For ETT-HYB use the data from 2000-2001 (Some problem with data in
        %2002
        k = nanmean(squeeze(OBSdH(2,:,:)-OBSdH(3,:,:)),1);
        k=circshift(k,[0,5]); % by circshifting the data 5 hoours and starting the center at 1 (not 0.5 hours) we achive UT to LT
        obs_index = 2;
        data1(:,1) = mjd;
        data1(:,2) = reshape(squeeze(OBSdH(2,:,:)-OBSdH(3,:,:))',[1,1096*24]);
        load c:\manoj\projects\eej\ettcm4.mat index sq_est;
        data1(:,1) = data1(:,1) + ((78-82.5)/15)/24; %NOTE this correction is for the Solar and IST difference between 78 (ETT) and 82.5 (Allahabad - IST)
        xtime = 1:24;
        obslon = 82.5; %Oservatory Longitude of ETT is given as 82.5 as the IST is referenced to this longitude
        st = 12.8; % The bins are still centered around ETT 
        en = 132.8;
        plotoptions = 'b*-';

    case 2,%ETT-ABG
        k = nanmean(squeeze(OBSdH(2,:,:)-OBSdH(1,:,:)),1);
        k=circshift(k,[0,5]); % by circshifting the data 5 hoours and starting the center at 1 (not 0.5 hours) we achive UT to LT
        obs_index = 2;
        data1(:,1) = mjd;
        data1(:,2) = reshape(squeeze(OBSdH(2,:,:)-OBSdH(1,:,:))',[1,1096*24]);
        data1(:,1) = data1(:,1) + ((78-82.5)/15)/24; %NOTE this correction is for the Solar and IST difference between 78 (ETT) and 82.5 (Allahabad - IST)
        xtime = 1:24;
        obslon = 82.5;
        st = 12.8;
        en = 132.8;
        plotoptions = 'b*-';

    case 3,%TIR-HYB
        k = nanmean(squeeze(OBSdH(4,:,:)-OBSdH(3,:,:)),1);
        k=circshift(k,[0,5]); % by circshifting the data 5 hoours and starting the center at 1 (not 0.5 hours) we achive UT to LT
        obs_index = 2;
        data1(:,1) = mjd;
        load c:\manoj\projects\eej\tircm4.mat index sq_est;
        data1(:,2) = reshape(squeeze(OBSdH(4,:,:)-OBSdH(3,:,:))',[1,1096*24]);
        data1(:,1) = data1(:,1) + ((77.85-82.5)/15)/24; %NOTE this correction is for the Solar and IST difference between 78.5 (TIR) and 82.5 (Allahabad - IST)
        xtime = [1:24];
        obslon = 82.5;
        st = 12.8;
        en = 132.8;
        plotoptions = 'b*-';

    case 4,%TIR-ABG
        k = nanmean(squeeze(OBSdH(4,:,:)-OBSdH(1,:,:)),1);
        k=circshift(k,[0,5]); % by circshifting the data 5 hoours and starting the center at 1 (not 0.5 hours) we achive UT to LT
        obs_index = 2; 
        data1(:,1) = mjd;
        load c:\manoj\projects\eej\tircm4.mat index sq_est;
        data1(:,2) = reshape(squeeze(OBSdH(4,:,:)-OBSdH(1,:,:))',[1,1096*24]);
        data1(:,1) = data1(:,1) + ((75-82.5)/15)/24; % This assumes 75 (between TIR and ABG)
        xtime = [1:24];
        obslon = 82.5;
        st = 12.8;
        en = 132.8;
        plotoptions = 'k.-';

    case 5,%HUA-FUQ
        k = nanmean(squeeze(OBSdH(5,:,:)-OBSdH(6,:,:)),1);
        k = circshift(k,[0,19]);
        xtime = 0.5:1:23.5;       
        obs_index = 2;
        data1(:,1) = mjd;
        data1(:,2) = reshape(squeeze(OBSdH(5,:,:)-OBSdH(6,:,:))',[1,1096*24]);
        data1(isnan(data1(:,2)),2) = 9999; 
        obslon = 285;
%        st = -142;
%        en = -8;
        st  = -139.5
        en = -10.5
        load c:\manoj\projects\eej\huacm4.mat index sq_est;
        plotoptions = 'k*-';
        
    case 6,%MBO-GUI
        k = nanmean(squeeze(OBSdH(7,:,:)-OBSdH(8,:,:)),1);
        k = circshift(k,[0,23]); % 23 ~ = 343/15
        %xtime = 0.5:1:23.5;       
        xtime = 1:24;
        data1(:,1) = mjd;
        data1(:,2) = reshape(squeeze(OBSdH(7,:,:)-OBSdH(8,:,:))',[1,1096*24]);
        data1(isnan(data1(:,2)),2) = 9999; 
        obs_index = 2;
        obslon = 343;
        st = -81.5;
        en = 47.5;
        plotoptions = 'g*-';
        load c:\manoj\projects\eej\mbocm4.mat index sq_est;
    case 7, %AAE-QSB
        k = nanmean(squeeze(OBSdH(15,:,:)-OBSdH(10,:,:)),1);
        k = circshift(k,[0,2]); % 2 ~ = 37.5/15
        xtime = 1:24;       
        obs_index = 2;   
        data1(:,1) = mjd;
        data1(:,2) = reshape(squeeze(OBSdH(15,:,:)-OBSdH(10,:,:))',[1,1096*24]);
        data1(isnan(data1(:,2)),2) = 9999; 
        obslon = 37.5;
        st = obslon - (incrlat*6+incrlat/2);
        en = obslon + (incrlat*6+incrlat/2);
        plotoptions = 'k*-';
        load c:\manoj\projects\eej\aaecm4.mat index sq_est;
    case 8,%GUA-CBI
        k = nanmean(squeeze(OBSdH(11,:,:)-OBSdH(12,:,:)),1);
        k = circshift(k,[0,9]); % 9 ~ = 144/15 (the remaining would be taken care by shifting time axis to 1:24 (instead 0.5:1:23.5)
        data1(:,1) = mjd;
        data1(:,2) = reshape(squeeze(OBSdH(11,:,:)-OBSdH(12,:,:))',[1,1096*24]);
        data1(isnan(data1(:,2)),2) = 9999; 
        xtime = 1:24;       
        obs_index = 2;  
        obslon = 143;
        st = obslon - (incrlat*6+incrlat/2);
        en = obslon + (incrlat*6+incrlat/2);
        plotoptions = 'r*-';
        load c:\manoj\projects\eej\guacm4.mat index sq_est;
    case 9,%YAP=BIK
        %FOR YAP, swith on robust processing -> few very very bad data
        %points !
        %YAP data are minute means 
%option 1 minute means
%         load c:\manoj\projects\eej\YAPMNBIK;
%         k = nanmean(YAPMNBIK);
%         k = circshift(k,[1,554]); %554 = (138.5/15)*60 minutes (UT shift in minutes)
%         xtime = (linspace(0,24,1440));       
%         data1(:,1)= yapdate;
%         data1(:,2) = yapmnbik;
% Option 2 hourly means
        load c:\manoj\projects\eej\YAPMNBIKHR;
        k = nanmean(YAPMNBIKHR);
        k = circshift(k,[1,9]); %9 = (138.5/15)hours (UT shift in minutes)
        xtime = 1:24;
        data1(:,1) = yaphrdate1;
        data1(:,2) = yapmnbikhr;
% COMMON STUFF
        obs_index = 2;
        obslon = 138.5;
        st = obslon - (incrlat*6+incrlat/2);
        en = obslon + (incrlat*6+incrlat/2);
        plotoptions = 'c*-';
        ncoeff = 9;
   case 10,%PND-HYB
        k = nanmean(squeeze(OBSdH(13,:,:)-OBSdH(3,:,:)),1);
        k=circshift(k,[0,5]); % by circshifting the data 5 hoours and starting the center at 1 (not 0.5 hours) we achive UT to LT
        obs_index = 2;
        data1(:,1) = mjd;
        data1(:,2) = reshape(squeeze(OBSdH(13,:,:)-OBSdH(3,:,:))',[1,1096*24]);
        %data1(:,1) = data1(:,1) + ((80-82.5)/15)/24; %NOTE this correction is for the Solar and IST difference between 80 (PND) and 82.5 (Allahabad - IST)
% Figure 13 requires that the above time correction not to be done ! (18.
% APR 2006)
        xtime = 1:24;
        obslon = 82.5; %79 degrees 55 minutes
%         st = obslon - (incrlat*6+incrlat/2);
%         en = obslon + (incrlat*6+incrlat/2);
        st = 15; % The bins are still centered around ETT 
        en = 135;
        load c:\manoj\projects\eej\pndindex.mat sq_est index
%        load c:\manoj\projects\eej\ettcm4.mat index sq_est;
        plotoptions = 'b*-';
   case 11,%ETT-PND
        k = nanmean(squeeze(OBSdH(2,:,:)-OBSdH(13,:,:)),1);
        k=circshift(k,[0,5]); % by circshifting the data 5 hoours and starting the center at 1 (not 0.5 hours) we achive UT to LT
        obs_index = 2;
        data1(:,1) = mjd;
        data1(:,2) = reshape(squeeze(OBSdH(2,:,:)-OBSdH(13,:,:))',[1,1096*24]);
        load c:\manoj\projects\eej\ettcm4.mat index sq_est;
        data1(:,1) = data1(:,1) + ((78-82.5)/15)/24; %NOTE this correction is for the Solar and IST difference between 78 (ETT) and 82.5 (Allahabad - IST)
        xtime = 1:24;
        obslon = 82.5; %Oservatory Longitude of ETT is given as 82.5 as the IST is referenced to this longitude
        st = 12.8; % The bins are still centered around ETT 
        en = 132.8;
        plotoptions = 'b*-';
   case 12,%TIR-PND
        k = nanmean(squeeze(OBSdH(4,:,:)-OBSdH(13,:,:)),1);
        k=circshift(k,[0,5]); % by circshifting the data 5 hoours and starting the center at 1 (not 0.5 hours) we achive UT to LT
        obs_index = 2; 
        data1(:,1) = mjd;
        load c:\manoj\projects\eej\tircm4.mat index sq_est;
        data1(:,2) = reshape(squeeze(OBSdH(4,:,:)-OBSdH(13,:,:))',[1,1096*24]);
        data1(:,1) = data1(:,1) + ((75-82.5)/15)/24; % This assumes 75 (between TIR and ABG)
        xtime = [1:24];
        obslon = 82.5;
        st = 12.8;
        en = 132.8;
        plotoptions = 'k.-';   
   case 13, %AAE-ELT
        k = nanmean(squeeze(OBSdH(15,:,:)-OBSdH(14,:,:)),1);
        k = circshift(k,[0,2]); % 2 ~ = 37.5/15
        xtime = 1:24;       
        obs_index = 2;   
        data1(:,1) = mjd;
        data1(:,2) = reshape(squeeze(OBSdH(15,:,:)-OBSdH(14,:,:))',[1,1096*24]);
        data1(isnan(data1(:,2)),2) = 9999; 
        obslon = 37.5;
        st = obslon - (incrlat*6+incrlat/2);
        en = obslon + (incrlat*6+incrlat/2);
        plotoptions = 'k*-';
        load c:\manoj\projects\eej\aaecm4.mat index sq_est;        
     otherwise,
        display('Error');
        break;
        return;

end;


%INFORMATION
%     The EEJ parameters in the file EEJ.TXT (version JUNE 2006)
% 1.  FDAY
% 2.  geographic longitude of the magnetic dip equator
% 3.  geographic latitude of the magnetic dip equator
% 4.  local_time at the magnetic dip equator
% 5.  max current (old definition, assuming that we can estimate the zero line)
% 6.  av. current difference to north and south troughs 
% 7.  current difference to south trough 
% 8.  current difference to north trough 
% 9.  magnetic field strength at 108 km altitude
% 10. Kp
% 11. F10.7
switch satdata,
    case 1,
        eej_index = 5;
    case 2,
        eej_index = 6;
    otherwise,
        display('Error');
        break;
        return;
end;        

if sqcor == 1,
    data1(index,obs_index) = data1(index,obs_index)-(sq_est(:,1)-sq_est(:,2));
   % data1(index,obs_index) = (sq_est(:,1)-sq_est(:,2));
end;

[p,s] = polyfit(xtime,k,ncoeff);%In case of minute data change the degree to 20
%----------------------------

ndd = zeros([1,25]);
nloop = 1;
lto = [];
ltc = [];
xdata = [];
nbindata = [];
ceff = [];
cerr = [];
nbin = 1;

for ki = st:incrlat:en,
L = data(:,2) > ki & data(:,2) < ki+incrlat & data(:,1) < maxdate  &  data(:,1) > mindate; %The last filter is to select only data before Dec 31, 2002
IndiaE = data(L,:);
ndata = 1;
ind1 = [];
ind2 = [];
ltt1 = [];
ltt2 = [];
nnodata=1;
for i = 1:length(IndiaE(:,1)),
    ind = findnearest(IndiaE(i,1),data1(:,1));
    [lt1,dummy1] = champ_lt(IndiaE(i,1),IndiaE(i,2),obslon);
    [dummy2,lt2] = champ_lt(data1(ind,1),IndiaE(i,2),obslon); % This is because observatory FDAY ~= CHAMP FDAY
	ltt1(i) = lt1;
    ltt2(i) = lt2;
    ltt3(i) = dummy1;
    ltt4(i) = dummy2;
    obst(i) = data1(ind,1);
    chpt(i) = IndiaE(i,1);
    kp1 = IndiaE(i,10);
    cmax = IndiaE(i,5);
    if lt2 >= obslt(1) & lt2 <= obslt(2) & lt1 >= obslt(1) & lt1 <= obslt(2) ...
            & abs(data1(ind,obs_index)) < 300 & data1(ind,1) < max(data(:,1)) ...
            & kp1 <= kp_thrashold & cmax >= 0.03, %The second last fliter take care the nodata (9999) points
            ind1(ndata) = ind;
            f_ch(ndata) = polyval(p,lt1); %Getting the expexted EEJ strength from poly coefficients for champ LT
            f_ob(ndata) = polyval(p,lt2); %Getting the expexted EEJ strength from poly coefficients for obs LT
%             f_ch(ndata) = interp1(1:24,k,lt1,'spline'); %Found spline interpolation working better 18.1.06 !
%             f_ob(ndata) = interp1(1:24,k,lt2,'spline');
%             fprintf('Lon = %5.1f i = %d f_ch = %f f_ob = %f lt1 = %f lt2 = %f kp = %3.1f\n', ki,i,f_ch(ndata),f_ob(ndata),lt1,lt2,data1(ind,end));
            ndata = ndata+1;
        else,
            ind2(nnodata) = i;
            nnodata = nnodata + 1;
        end;
end;

if length(ind2) > 0,
    IndiaE(ind2,:) = [];
    obst(ind2) = [];
    chpt(ind2) = [];
end;

if length(ind1) > 1,
    obsE = data1(ind1,:);
    X = IndiaE(:,eej_index);
    if ltcor == 1,
        Y = obsE(:,obs_index).*(f_ch./f_ob)';
        elseif ltcor == 0,
        Y = obsE(:,obs_index)';
    end;
    
    if robft == 1,    
%    Robust outlier rejection  -- START
     if length(X) > 3,
     [rs,rstat] = robustfit(X,Y);
     residuals = abs(rstat.resid);
     LLL = residuals == max(residuals); %Just remove the most noisy point
     X(LLL) = [];
     Y(LLL) = [];
     end;
    % Robust outlier rejection  --- END
    end;

    [cof,err] = corrcoef(X,Y);  % correlation & Correction for observatory data 
    ceff(nloop) = cof(2,1);
    cerr(nloop) = err(2,1);
    fprintf('%4.3f ',cof(2,1));
    xdata(nloop) = ki+incrlat/2;
    ndd(nloop) = ndata-1;
    lto(nloop) = min(ltt2);
    ltc(nloop) = max(ltt2);
    nloop = nloop +1;

    %----------- PLOT SCATTER
%subplot(2,3,nloop-1);
if x1data(nbin) == 0,

figure(2);
% hold on;
plot(X,Y,[plotcol 's'],'MarkerFaceColor',plotcol);
set(gca,'FontSize',16);
grid on
hold on;
% %axis([0 0.5 0 200]);
h=xlabel('A/m');
% % %h=xlabel(['Longitude Bin ' sprintf('%4.1f ',ki) '^\circ  '  sprintf('%4.1f ',ki+incrlat) '^\circ']);
set(h,'FontSize',16);
ylabel('\DeltaH');
[c,sfit] = polyfit(X,Y,1);
m = polyval(c,[min(X),max(X)]);
plot([min(X),max(X)],m,plotopt,'LineWidth', 2);
text(0,25,sprintf('dH = %5.1f * I + %5.2f',c(1),c(2)),'FontSize',16,'Color',plotcol);
end;


%-----------------END
else,
    x1data(nbin) = NaN;
     
end;
nbin = nbin + 1;
clear ind f_ch f_ob ind1 ind2 ndata nnodata IndiaE obsE ltt1 ltt2 ltt3 ltt4 obst chpt X Y;
end;

x1data(isnan(x1data)) = [];
fprintf('\n');
figure(1);
hold on;
errorbar(x1data,ceff,cerr/2,plotopt,'LineWidth',2);
%axis([-80 80 -.6 1.4]);
grid on;
hold on;